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Surface and ferroelectric properties of PbTiOs thin films are investigated using an interatomic 
potential approach with parameters computed from first-principles calculations. We show that a 
model developed for the bulk describes properly the surface properties of PbTiOs. In particular, 
the antiferrodistortive surface reconstruction, recently observed from X-ray scattering, is correctly 
reproduced as a result of the change in the balance of long-range Coulombic and short-range inter- 
actions at the surface. The effects of the surface reconstruction on the ferroelectric properties of 
ultrathin films are investigated. Under the imposed open-circuit electrical boundary conditions, the 
model gives a critical thickness for ferroelectricity of 4 unit cells. The surface layer, which forms the 
antiferrodistortive reconstruction, participates in the ferroelectricity. A decrease in the tetragonality 
of the films leads to the stabilization of a phase with non- vanishing in-plane polarization. A peculiar 
effect of the surface reconstruction on the in-plane polarization profile is found. 
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I. INTRODUCTION 

One fascinating feature of perovskites is that they ex- 
hibit a large variety of structural phase transitions. The 
variety of transition behavior and low-temperature dis- 
torted structures depend on the individual compound. 
Among the perovskites one finds ferroelectric (FE) crys- 
tals such as BaTiC>3, KNbC>3 (displaying three phase 
transitions), and PbTiC>3 (displaying only one transi- 
tion), antiferroelectrics such as PbZrC>3, and materials 
such as SrTiOs that exhibit other non-polar antiferrodis- 
tortive (AFD) instability involving the rotation of the 
oxygen octahedra 0. 

With the rapidly advancing miniaturization of ferro- 
electric devices and the use of thin films, attention is 
focusing on the role played by surfaces and interfaces in 
the overall performance of the materials. The effects of 
surfaces on the structural phase transitions have become 
an issue of significant importance because the surface can 
affect the structural behavior of the perovskites by modi- 
fying the strength of various instabilities. Lead titanate is 
a clear example. In PbTiC>3 (PT), ferroelectricity is due 
to the condensation of a Tis unstable phonon in which 
the oxygen octahedra shift against the Pb sub-lattice. 
The ground state consists of shifts along (001) accompa- 
nied by a tetragonal lattice strain which stabilizes this 
direction. However, complete phonon dispersions for 
the ideal perovskite structure PbTiOs as determined by 
density functional calculations show ferroelectric (r^) 
and also rotational antiferrodistotive (R25 type) instabil- 
ities. Q In bulk, the FE and AFD instabilities compete 
with each other and the FE lattice distortion suppresses 
the AFD distortion, but the proximity to a surface mod- 
ifies that balance. Recently, an antiferrodistortive recon- 
struction of the PbTiOs (001) surface has been found 
using grazing incidence X-ray scattering. The atomic 



structure of the surface consist of a single layer of an 
AFD structure with oxygen cages rotated by 10° around 
the [001] axis through the Ti ions. Latter on, the AFD 
reconstruction at the PbO-terminated (001) surface was 
reproduced by Bungaro and Rabe using ab-initio calcu- 
lations. They also found for in-plane polarized films 
that the FE and AFD distortions coexist in the proxim- 
ity of the surface. Regarding device applications, out-of- 
plane polarized films are more relevant. Therefore, it is 
worth investigating also the effects of the surface recon- 
struction on the out-of-plane ferroelectricity, since it was 
speculated that the AFD reconstruction may be related 
to the dead layer often postulated to explain the observed 
decrease in ferroelectric character in thin films. Q 

The realistic simulation of PT thin- films is a theoretical 
challenge due to the interplay of polar and non-polar in- 
stabilities at the surface. Although first-principles meth- 
ods are extremely precise and include thoroughly elec- 
tronic effects, they are quite computer demanding. So, 
these calculations are restricted to the investigation of 
zero-temperature properties of perovskites involving a 
rather small number of atoms. For the study of the 
thermal behavior, or to handle larger system sizes, other 
methods are necessary. During the last years, the ef- 
fective Hamiltonian approach has been used for the in- 
vestigation of thin film properties H,SH>I3- However, 
for the simulation of Pb-based perovskite thin-films the 
coexistence of rotational and FE distortions should be 
explicitly considered. 01 Atomistic simulations based on 
interatomic potentials can account naturally for the pres- 
ence of competing instabilities. However, the validity 
of any atomistic simulation study depends to consider- 
able extent on the quality of the potentials used. The 
interatomic potential approach firmly grounded by hav- 
ing its parameters fitted against results of first-principles 
calculations is a promising technique for investigating 
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thin-film properties of Pb-based perovskites. In this pa- 
per we investigate the ferroelectric properties of PT thin 
films using a shell model developed completely from first- 
principles calculations. We investigate first if the model 
developed for the bulk describes properly surface prop- 
erties. Then the ferroelectric behavior of ultrathin films 
under open-circuit electrical boundary conditions are in- 
vestigated. 



II. MODEL AND COMPUTATIONAL DETAILS 

The approach we follow in this work is based on the 
atomistic modeling using interatomic potentials. We in- 
deed found that the shell model approach does provide 
computationally efficient and confident methodology for 
the simulation of ferroelectric perovskites, including bulk 
properties of pure crystals 0, 0, 0] and solid solu- 
tions EL superlattices ^| , surface and thin films prop- 
erties. In the shell model, each atom is represented 
by a massive core coupled to a massless shell, and the 
relative core-shell displacement describes the atomic po- 
larization. The model contains coulombic long-range and 
pairwise short-range interactions. In the present study 
we use the model developed in Ref |13| which contains 
4 th order core-shell couplings and short-range interac- 
tions described by two different types of potentials. A 
Rydberg potential (V(r) = (a + br) exp(— r / p)) is used 
for the Pb-Ti, Pb-0 and Ti-0 pairs, and a Buckingham 
potential (V(r) = aexp(— r/p) + cr~ 6 ) is used for 0-0 
interactions. The model was developed completely from 
first-principles calculations within the Local Density Ap- 
proximation (LDA). The parameters were adjusted to 
reproduce a wide number of LDA results, including lat- 
tice dynamics and total energy calculations. The result- 
ing model was able to reproduce delicate properties of 
PT in good agreement with LDA results. However, as 
a consequence of the LDA, volume and volume depen- 
dent properties are underestimated with respect to ex- 
perimental values. For example, the model gives a tetrag- 
onal ground state with a lattice parameter a=3.859A, a 
tetragonal distortion c/a= 1.043, and a spontaneous po- 
larization P=54/xC/cm 2 , while the experimental data are 
a=3.90l, c/a=1.065, and P=75/iC/cm 2 . The underes- 
timation of the static structural properties is translated 
via the adjusted model to the finite temperature behav- 
ior. Molecular dynamics simulations showed a cubic- 
tetragonal transition at Tc — 450K, 300K below the 
experimental value. Nevertheless, the qualitative tem- 
perature behavior of lattice parameters and polarization 
was correctly reproduced. See Reference [l3j for more 
details. 

The investigation of the surface properties is carried 
out using an isolated slab geometry with periodic bound- 
ary conditions in the x-y plane. The long ranged electro- 
static energy and forces are calculated by a direct sum 
method |l7| . The equilibrated zero-temperature struc- 
ture of the slabs was determined by a zero-temperature 



quench until the force on each individual ion was less than 
0.001 eV / A. To mimic the two dimensional clamping and 
straining of the film due to the presence of a substrate, 
we force the simulation cell to be square in the x-y plane. 



III. RESULTS AND DISCUSSION 

A. Surface effects in non-polar slabs. Validation of 
the model 

The parameters of the model were derived to describe 
perfect crystal properties, and it is not guaranteed that 
the same potentials are suitable to describe the surface 
properties of PT. In the bulk, the underlying potential 
surface is a result of a delicate balance between long range 
Coulomb, short range, and core-shell interactions. Such 
balance could change significantly at the surface, where 
the atomic environment changes due to the discontinuity 
of the crystal. So, it is important to check if the model 
developed for the bulk proves also successful for describ- 
ing surface properties. 

A first-principles study for (001) surfaces of cubic PT 
has been performed by Meyer, Padilla and Vanderbilt 
using c(l x 1) surface periodicity [18|] As a first step, 
we take these first-principles calculations as benchmark 
results to compare with. To this end we determined the 
atomic equilibrium position for both PbO- and TiO- sym- 
metrically terminated 7-layer periodic slabs. We set the 
in-plane lattice parameter to the cubic equilibrium value 
yielded by the model in bulk, a=3.887A, and we relaxed 
the atomic positions from the ideal positions only in the 
direction perpendicular to the surface. The results ob- 
tained for the atomic relaxations for both slabs are listed 
in Table ^ whereas the values for the interlayer dis- 
tances and rumplings are listed in Table [HI Ab-initio 
results are shown in parenthesis as reference. The model 
reproduces satisfactorily the relaxation direction of the 
atoms, the change in the interlayer distances, and atomic 
rumplings. For the case of the PbO surface, the equilib- 
rium termination for PT |18| . the numerical values are 
in very good agreement with the ab-initio ones. The av- 
erage surface relaxation energy obtained by the model 
AE re i ax = Q.237eV is also in quite good agreement with 
the ab-initio result of 0.210 eV. 

The recent pseudopotential study of the AFD surface 
reconstruction in PbTiOa slabs Jf| provides reliable ad- 
ditional information to validate the model. In Table lllll 
we show structural parameters of the relaxed c(2 x 2) 
11-layer PbO terminated slab. The overall agreement 
with the ab-initio results is very good. The model re- 
produces the main differences between the c(2 x 2) and 
(1 x 1) surface periodicities, that is the presence of an 
AFD reconstruction with a reduction of rumplings and 
interlayer distances. Regarding the AFD surface recon- 
struction, the rotation angles are slightly overestimated 
with respect to the pseudopotential calculations; the ex- 
perimentally determined rotation angle for the surface 



3 



Atom 


5 Z 


Atom 


S z 


Pb(l) 


-3.45 (-4.36) 


Ti(l) 


-4.17 (-3.40) 


O ra (l) 


-0.27 (-0.46) 


0/(1) 


-2.94 (-0.34) 


Ti(2) 


3.19 ( 2.39) 


0//(l) 


-2.94 (-0.34) 


0/(2) 


1.65 ( 1.21) 


Pb(2) 


1.58 ( 4.53) 


0// (2) 


1.65 ( 1.21) 


Om(2) 


-0.34 ( 0.43) 


Pb(3) 


-0.78 (-1.37) 


Ti(3) 


-0.98 (-0.92) 


0///(3) 


0.33 (-0.20) 


0/(3) 


0.56 (-0.27) 






0//(3) 


0.56 (-0.27) 



TABLE I: Atomic relaxations perpendicular to the surface 
(8z) of the PbO (left panel) and the Ti02 (right panel) termi- 
nated surface in the cubic phase. S z are given as percent of the 
theoretical unit cell parameter a. For comparison, ab-initio 
results (Ref are shown in parentheses. 



PbO surface 


TiO surface 




-4.04 (-4.2) 


-4.17 (-4.4) 


Ad 23 


2.34 ( 2.6) 


+0.83 (+3.1) 


Ad 34 


-0.02 (-0.8) 


-0.21 (-0.6) 


V\ 


3.02 (3.9) 


1.22 (3.1) 


Vl 


1.85 (1.2) 


1.93 (4.1) 


V3 


0.96 (1.2) 


1.54 (0.7) 



TABLE II: Change in the interlayer distance (Ad) and layer 
rumpling (v) for the PbO (left panel) and the Ti02 (right 
panel) terminated surface in the cubic phase, given as percent 
of the theoretical unit cell parameter a. For comparison, ab- 
initio results (Ref ^3) are shown in parentheses. 



layer is 10° The energy gain associated with the oc- 
tahedra rotations is 0.32 eV per surface unit cell, which 
is approximately one order of magnitude larger than the 
bulk ferroelectric well depth. 

It was pointed out that the enhancement of the AFD 
distortion at the PbO surface is a consequence of the 
particular chemistry of Pb, due to its tendency to move 
off-center and form strong covalent Pb-0 bonds We 
have demonstrated that this electronic effect, which has 
been showed to be an important factor for ferroelectricity 
in PbTiC>3 2], can be mimicked at the atomic level by 
a shell model. In this less-sophisticated approach, the 
AFD surface reconstruction is obtained just as a result 
of the change in the balance of long-range Coulombic and 
short-range interactions at the surface. 



Interlayer 


Rumpling 


Angle 


Adl2 


-3.3 (-3.4) 


V\ 


-1.4 (-1.4) 


02 


Ad 23 


+2.1 (+2.9) 


V-2 


+1.3 (+0.9) 


04 


Ad 34 


-0.1 (-0.9) 


Vi 


-1.4 (-2.0) 


06 


Ad 45 


+0.3 (+0.4) 


U A 


+0.5 (+0.4) 


Qbulk 


Ad 56 


+0.1 (-0.1) 


l>5 


+0.2 (-0.2) 





B. Ferroelectric behavior in polar slabs 

We use the atomistic model to explore zero- 
temperature ferroelectric properties in PT ultrathin 
films. We consider only PbO-terminated surfaces be- 
cause they are the most stable ones. 0] The simulation 
cell contains 10 x 10 unit cells with periodic boundary 
conditions along the x — y plane. The in-plane lattice pa- 
rameter was set and clamped to the model equilibrium 
value for the tetragonal bulk a=3.859A Since this value 
is close to the lattice parameter of SrTiC>3, this clamped 
condition simulates the strain effect of a SrTiC>3 sub- 
strate. 

We have performed standard-atomic relaxation meth- 
ods to determine the zero temperature structure in slabs 
from 2 to 10 unit cells thick that contains from 1200 to 
5200 atoms, respectively. In each case, the relaxation 
started with atoms slightly displaced from their ideal 
cubic positions. After a first quenching, the system is 
iteratively warmed-up and quenched again in order to 
eliminate possible metastable configurations. To analyze 
the results, we define the cell parameter c as the distance 
between two consecutive PbO planes, and the local po- 
larization is defined as the polarization of a Ti-centered 
unit cell. 

A decrease in tetragonality is an important effect pro- 
duced by the presence of the surface. Figure 1 provides 
detailed microscopic information about the average cell- 
by-cell tetragonal distortion c/a for a 10-unit cells thick 
slab. It can be seen that the tetragonality of all unit 
cells in the slab is considerable lower than in the bulk 
(c/a = 1.043 for the bulk). The stronger reduction is 
observed for the surface cells which are practically cubic 
(c/a = 1.005) while the tetragonality increases towards 
the interior of the slab reaching the value c/a — 1.032. As 
the imposed in-plane lattice parameter corresponds with 
the one of ferroelectric bulk PT, we can expect that the 
tetragonality gradually converges to the bulk value when 
the film thickness increases. In fact, we have obtained 
that, while the tetragonality of the outermost layer prac- 
tically does not depend on film thickness (surface relax- 
ation effects are predominant for these layers), the inner 
cells display a thickness dependence. Figure 2 shows the 
evolution of the c/a ratio with film thickness for cells in 
the middle of the slab. We can roughly estimate from the 
extrapolation of Figure 2 that d = %QA (s=s 22 unit cells) is 
_the lower-bond estimate of the slab thickness requiered 
io reach the tetragonality of the bulk in the center of 



13.8 (H.4) the 

film. The strong thickness dependence of tetragonal 



-1.7 (-2.9). 
8.0 (3.9) 
6.0 (3.3) 



TABLE III: Change in the interlayer distance (Ad), layer 
rumpling (y) and rotation angle of oxygen octahedra (0) for 
the PbO terminated surface at a\, u ik in a 11-layer non-polar 
PbTiOs slab. For comparison, ab-initio results (Ref. [f|) are 
shown in parentheses. 



ity is in qualitative agreement with recent experimental 
measurements which showed that the average c/a ratio 
decreases substantially for films thinner than 200A . 

In slabs under ideal open-circuit electrical boundary 
conditions, as we are considering here, out-of-plane fer- 
roelectricity manifests above a critical thickness, through 
the formation of stripe domains of alternating polarity. 
These domains are effective in neutralizing the depolar- 
ized field and stabilizing the ferroelectric phase, even 
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FIG. 1: Average cell-by-cell tetragonal distortion c/a for a 
10-unit cells thick slab. 
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FIG. 2: Thickness dependence of the tetragonal distortion (for 
the inner cells) and the spontaneous out-of-plane polarization 
of the ferroelectric domains. 



without electrodes. 0, The observation of stripe 
domains below the ferroelectric transition in epitaxial 
PbTiC>3 ultrathin films on SrTiC>3 substrates was re- 
ported from X-ray scattering |20j |. and the ferroelectric 
phase was found to be stable for thickness down to 3 
unit cells |2l|. Figure 3 shows the cell- by-cell out-of- 
plane polarization profile of a chain perpendicular to the 
surface of slabs of 4 and 5 unit cell width. It is clear that 
under the imposed stress and electrical boundary condi- 
tions, the model gives a critical thickness (d c ) for ferro- 



electricity of 4 unit cells (d c is defined as the maximun 
thickness at which polarization is zero). In the 4-cell 
slab the surface cells develop a small inward polarization 
of ps 3fiC/cm 2 due to surface atomic relaxations (see 
Table UJ) while the polarization of the two inner cells is 
practically zero. In the 5-cell slab, individual chains have 
a net out-of-plane polarization but the average polariza- 
tion of the slab is zero due to the development of stripe 
domains (the evolution of the polarization at the center 
of the slab with film thickness is shown in Figure 2). The 
inner cells display a local polarization of w 18/iC/cm 2 , 
that is 0.33 of the bulk value. The surface cells present an 
even stronger reduction of polarization. It is very impor- 
tant to point out, however, that the polarization of the 
surface cells in the 5-cell ferroelectric slab is considerable 
larger than in the 4-cell non-ferroelectric one, it increases 
narsksfrom ~ 3/zC/cto 2 to w 10/iC/cm 2 . This polarization en- 
hancement indicates that the surface layer, which forms 
an AFD reconstruction, also participates in the ferroelec- 
tricity. We note that the value of d c obtained with the 
model is two unit cells larger than the experimental one 
reported in Ref. |2l|. This difference could arise from 
the fact the model was developed to reproduce the LDA 
ground state of bulk PT, and consecuently the ferroelec- 
tric instability is underestimated. Another possibility is 
that we are not explicity simulating the presence of a 
substrate (although we are considering its strain effects) 
which modifies the boundary condition at one face of the 
film. Atomistic investigations of the substrate effects on 
the critical thickness are in progress. 

In the simulations, the in-plane lattice parameter was 
set to the tetragonal bulk value a=3.859A, so a vanish- 
ing in-plane polarization component could be expected. 
However, this is not the case and the slabs develop a 
polarization component along the (110) direction. For 
the bulk, the coupling of the ferroelectric soft mode 
with strain is the main driving force for the stabilization 
of the (OOl)-polarized tetragonal phase over the (111)- 
polarized rhombohedral one Q. In films, the decrease 
in tetragonality favors the stabilization of a phase with 
a non- vanishing in-plane polarization component, above 
and also below the critical thickness. Figure 4 shows 
the change of the average in-plane polarization with film 
thickness. The reduction of tetragonality produces that 
the in-plane polarization increases when the film thick- 
ness decreases. Locally, a strong increment of the in- 
plane polarization at the surface is observed (see the pro- 
file in the inset of Figure 4) . The polarization decreases 
towards the interior of the slab. For thick enough films, 
it is expected that the interior of the slab reachs the van- 
ishing in-plane polarization of the bulk. 

We showed above that non-polar slabs displayed an 
AFD surface reconstruction. Polar slabs show a surface 
reconstruction with a similar rotation angle which is in- 
dependent on the slab thickness. However, one difference 
is that neighboring planes of octahedra along the rotation 
axis rotates in-phase instead of out-of-phase. In fact, for 
the polar slabs, 8 2 = 13.5°, 4 = 4.5° and 6 6 = 7.3° for 
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FIG. 3: Cell-by-cell out-of-plane polarization profile of a ran- 
domly chosen chain perpendicular to the surface for a slab of 
4 and 5 unit cell width. 



In summary, we have shown that the delicate balance 
of lattice strain and polar and non-polar instabilities, 
which is responsible for the observed tetragonal ground- 
state of bulk PbTiC>3 and the AFD structure of the sur- 
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Figure 3 -Sepliarsky 



the first, second and third layers, respectively. We note 
that unstable modes at the M and R points, which are 
associated with oxygen rotation instabilities, have simi- 
lar frequencies and compete to each other in the cubic 
bulk • For the bulk, the model yields to an energy dif- 
ference of only 2 meV per unit cell between the M and R 
distortions (A E = 1.2 meV per unit cell from ab-initio 
calculations |j|). The presence of the polarized surface 
affects these two competing structural instabilities stabi- 
lizing M-type distortions at the surface. 

Finally, we show a peculiar effect of the AFD surface 
reconstruction on the in-plane polarization profile. Al- 
though the average in-plane polarization is oriented along 
the (110) direction, unit cells are not uniformly polarized 
along that direction. This can be visualized by the top 
view of the surface showed in Figure 5. We will consider 
the displacements of the Pb atoms to make the descrip- 
tion simpler. It is clear from the figure that the Pb atoms 
displace mainly along (100) directions (see the arrows) 
in such a way that half of the Pb atoms are displaced 
along a (100) direction, and the other half along a (010) 
direction. Although the effect is much stronger at the 
surface, this local behavior is observed throughout the 
ultrathin film. This polarization profile, which produces 
an average in-plane polarization along the (110) direc- 
tion, is a consequence of the oxygen octahedra rotation 
plus the tendency of Pb to move off-center shortening 
Pb-0 bonds. 



FIG. 4: Evolution of the average in-plane polarization with 
film thickness. Inset: cell-by-cell in-plane polarization profile 
for a randomly chosen chain perpendicular to the surface in 
a slab of 10 unit cell width. 



face, is accurately simulated by an interatomic potential 
approach with parameters computed from first-principles 
calculations. In ultrathin films, a remarkable decrease in 
tetragonality leads to the stabilization of a phase with 
non-vanishing in-plane polarization, above and also be- 
low the critical thickness (d c ) for ferroelectricity. Under 
the imposed stress and electrical boundary conditions, 
d c was found to be 4 unit cells. The surface layer, which 
forms an AFD reconstruction, also participates in the 
ferroelectricity. 
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